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' Abstract: A chain regularization method is combined with special pur- 

pose computer hardware to study the evolution of massive black hole 
binaries at the centers of galaxies. Preliminary results with up to A = 
0.26 x 10 6 particles are presented. The decay rate of the binary is shown 
fN| . to decrease with increasing N, as expected on the basis of theoretical 

arguments. The eccentricity of the binary remains small. 
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Coalescence of binary supermassive black holes is potentially the strongest 
source of gravitational waves in the universe [1]. The coalescence rate is limited 
by the efficiency with which massive binaries can interact with stars and gas in a 
galaxy and reach the relativistic regime at separations of ~ 10 -3 pc. Exchange 
of energy between a binary black hole and stars should also leave observable 
traces in the stellar distribution, perhaps allowing us to infer something about 
the merger history of galaxies from their nuclear structure [2] . Henry Kandrup 
worked on this problem shortly before his death. In "Supermassive Black Hole 
Binaries as Galactic Blenders" [3], Kandrup et al. investigated the effects of 
a massive binary on the stellar orbits near the center of spherical and nearly 
spherical galaxies. They showed that the periodically- varying potential due to 
the binary, coupled with the fixed potential from the galaxy, was effective at 
inducing chaos in the stellar orbits, leading to diffusion in both energy and 
configuration space and to ejection of stars from the nucleus. This study was 
a complement to earlier studies based on scattering experiments [4,5] in which 
the potential of the galaxy was ignored. 

Another approach to the binary black hole problem is via direct A-body 
techniques [6-8]. This approach is computationally challenging because of the 
need to handle close interactions between the star- and black hole particles with 
high precision. In addition, large particle numbers are required to avoid the 
effects of spurious relaxation [9,10]. Here, we present preliminary results of N- 
body integrations of the binary black hole problem, in which close interactions 
between the black holes and stars are handled via the Mikkola- Aarseth chain- 
regularization algorithm [11,12]. Recently Aarseth [13] described an application 
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of a time transformed leapfrog scheme [14] to this problem. We prefer the 
chain algorithm since it has proved itself in numerous applications, including one 
very similar to the current problem [15]. We incorporate the chain algorithm 
into a general-purpose iV-body code by including the effects of nearby stars as 
perturbers to the chain. We present some preliminary results of binary black hole 
evolution computed via this algorithm on a special-purpose GRAPE-6 computer 
with particle numbers up to 0.26 x 10 6 . 

Our basic iV-body algorithm is an adaptation of the NB0DY1 code of Aarseth 
[16] to the GRAPE-6 special purpose hardware. The code uses a fourth-order 
Hermite integration scheme with individual, adaptive, block time steps [17]. For 
the majority of the particles, the forces and force derivatives were calculated via 
a direct-summation scheme using the GRAPE-6. 

Close encounters between the massive particles ("black holes"), or between 
black holes and stars, create prohibitively small time steps in such a scheme. To 
avoid this situation, we regularized the critical interactions as follows. Let r,, 
i = l, N be the position vectors of the particles. We first identify the subset 
of n particles to be included in the chain; the precise criterion for inclusion is 
presented below, but in the late stages of evolution, the chain always included the 
two black holes as its lowest members. We then search for the particle which is 
closest to either end of the chain and add it; this operation is repeated recursively 
until all n particles are included. Define the separation vectors R^ = r^+i — 
where r i+ i and r, are the coordinates of the two particles making up the ith 
link of the chain. The canonical momenta W ; corresponding to the coordinates 
Ri are given in terms of the old momenta via the generating function 

n-l 

5=^W r (r I+1 -4 (1) 

i=l 

Next, we apply KS regularization [18] to the chain vectors, regularizing only the 
interactions between neighboring particles in the chain. Let Qi and Pi be the 
KS transformed R^ and W; coordinates. After applying the time transormation 
5t = g5s, g = 1/L, where L is the Lagrangian of the system (L = T — U, where 
T is the kinetic and U is the potential energy of the system). We obtain the 
regularized Hamiltonian T = g(H(Qi,Pi) — E ), where E is the total energy of 
the system. The equations of motion are then 

, dT dT 

where primes denote differentiation with respect to the time coordinate s. Be- 
cause of the use of regularized coordinates, these equations do not suffer from 
singularities, as long as care is taken in the construction of the chain. 

Since it is impractical to include all iV particles in the chain, we must 
consider the effects of external forces on the chain members. Let Fj be the 
perturbing acceleration acting on the jth body of mass rrij. The perturbed 
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system can be written in Hamiltonian form by simply adding the perturbing 
potential: 



Only one chain was defined at any given time. At the start of the iV-body 
integrations, there was no regularization, and all particles were advanced using 
the variable-time-step Hermite scheme. The chain was "turned on" at the time 
when one of the particles (including possibly one of the black holes) achieved a 
time step shorter than t c hmin and reached a distance from one of the black holes 
smaller than r c h min . Each star inside r c h m in radius was then added to the chain, 
and the two black holes were always included. The values of t c h min and r c h min 
were determined by carrying out test runs; we adopted t c hmin ~ 10" 5 - 10" 6 
and r c hmin ~ 10~ 4 — 10~ 3 in standard iV-body units. 

The chain's center of mass was a pseudoparticle as seen by the ./V-body 
code and was advanced by the Hermite scheme in the same way as an ordinary 
particle. However, when integrating the trajectories of stars near to the chain, 
it is essential to resolve the inner structure of the chain. Thus for stars inside 
a critical r cr i t \ radius around the chain, the forces from the individual chain 
members were taken into account. The value of r crit i was set by the size of the 
chain to be r cr iti = XR c h with R c h the spatial size of the chain and A = 100. 

In addition, the equations of motion of the chain particles must include 
the forces exerted by a set of external perturber stars. Whether or not a 
given star was listed as a perturber was determined by a tidal criterion: r < 
Rcrit2 = (■m/mchain) 1/3 l m 1 inRch where m chain represents the mass of the chain, 
m is the mass of the star, and "y m in was chosen to be 10~ 6 ; thus r crit2 ~ 



The membership of the chain changed under the evolution of the system. 
Stars were captured into the chain if their orbits approached the binary closer 
than R c h . Stars were emitted from the chain if they got further from both of the 
black holes than 1.5R c h- The difference between the emission and absorption 
distances was chosen to avoid a too-frequent variation of the chain member- 
ship. When the last particle left the chain, the chain was eliminated and the 
integration turned back to the Hermite scheme, until a new chain was created. 

In our numerical experiments, the typical number of chain members was 
5 — 10. The number of perturber stars was typically 500 — 1000, and the number 
of stars inside r cr i t \ was 2000—5000. Using a minimum tolerance of 10 -12 for the 
chain's Bulirsch-Stoer integrator allowed us to reach typical relative accuracies 
of 10 -9 in the chain integration. The relative accuracy in the conservation of 
the total energy of the system was determined by the Hermite scheme. For all 
of our numerical simulation the relative error in the total energy was less than 
10~ 4 during the course of the integration. 
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For this set of experiments, we adopted Dehnen's [19] density law, 



p(r) 



(3 - j) M 
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with 7 = 0.5, to describe the initial stellar distribution. The initial positions 
and velocities of the N stars were generated from the steady state phase space 
density f(E) that reproduces the density law (4). Henceforth we adopt units 
such that the gravitational constant G, the total stellar mass M, and the Dehnen 
scale length a are equal to one. To this model we added two black hole particles 
each of mass 0.005. The black holes were placed symmetrically about the center 
of the galaxy, offset at a distance 0.1 from the center. The tangential velocities 
were set to ±0.16 yielding nearly circular initial orbits for the two black holes 
about the center of the galaxy. 

We integrated the above model with three different particle numbers: N = 
16384, 65536, and 262144; the latter is close to the maximum number of particles 
that can be handled in the GRAPE-6 memory. The integrations were carried 
out for 170 time units. Elapsed times were (2.5, 26, 105) hours for the three runs. 
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Figure 1. Evolution of the binary semi-major axis. 

The orbits of the two black holes initially decay, and at a time of roughly 30 
they form a bound pair. After this, the semimajor axis a of the binary shrinks as 
the two black holes interact with stars and eject them from the nucleus via the 
gravitational slingshot. The instantaneous decay rate is given approximately by 
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where p and a are the density and velocity dispersion of the stars, and H is 
a dimensionless constant of order 16 [4,5]. Because the stellar density changes 
with time as the binary ejects stars, the binary's decay rate can in general be a 
complicated function of time. Two limiting cases are of interest [9]. When the 
particle number N is small, gravitational encounters are able to scatter stars into 
the "loss cone" around the binary at a higher rate than they are ejected. The 
density near the binary remains approximately constant and the decay follows 
a -1 ~ t. This is the "full loss cone" regime. When N is large, encounters 
between stars are weak, and the binary's loss cone remains nearly empty. Decay 
of the binary is limited by the rate at which stars diffuse into the loss cone; 
since the diffusion time scales approximately as N, the binary decay follows 
a -1 ~ t/N. This is the "diffusion" regime. Real galaxies are expected to be in 
the diffusion regime [9]. 

Figure 1 shows the time evolution of the semi- major axis of the binary in 
each of the three integrations. The binary's energy evolves as an approximately 
linear function of the time but with a prefactor that depends on N; the approx- 
imate dependence is 

1 160t 

This is intermediate between the a -1 oc t dependence of the full loss cone limit, 
and the a -1 oc t/N dependence in the diffusion limit. We conclude that, for 
the particle numbers considered here, replenishment of the binary's loss cone is 
taking place but at a lower rate than the rate at which the loss cone is being 
emptied. Apparently, particle numbers in excess of ~ 10 6 are required if N- 
body integrations are to be completely in the diffusive regime characteristic of 
real galaxies. Such large particle numbers can in principle be handled with 
direct-summation codes like ours if coupled with parallel hardware [20]. 

Figure 2 shows the eccentricity evolution of the binary. The value of e at 
the time when the hard binary first forms, t ~ 30, is substantially different in 
the different integrations due presumably to finite-TV effects during the initial 
inspiral of the two black holes. Thereafter the eccentricity fluctuates in the case 
of the two small- N integrations, but decreases with time in the integration with 
the largest N. We compared the eccentricity evolution with the predictions of 
scattering theory. The rate of change of e is commonly written 

^ = ^lna" 1 (7) 
dt dt v ' 

where K = K(e,a) [4,5]. The functional form of K(e,a) is not well known; we 
adopted the expression given in [5] for a hard binary. We then wrote equation 
(7) as 

ei = ei-i + K (ej_i,aj_i)ln (^^j ( 8 ) 

where the subscript denotes the time step. Combining equation (8) with the 
./V-body results for a(t), Figure 1, we could then predict the expected evolution 
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Figure 2. Evolution of the binary orbital eccentricity. 



in e. The result for the largest-TV integration is shown as the dashed line in 
Figure 2. There is reasonable agreement, but the eccentricity evolution even in 
the largest-iV integration still exhibits substantial fluctuations. In this regard 
too, we are not yet in a regime where the evolution is similar to what would be 
expected in real galaxies. 

A recent iV-body study [13] found a much greater degree of eccentricity 
evolution, although the initial orbit of the binary was highly non-circular. 

As the binary decays, it ejects stars from the nucleus and lowers the nuclear 
density. Figure 3 shows initial and final density profiles for the three integrations. 
The net effect of a black hole binary on the stellar distribution is commonly 
measured in terms of the "mass deficit," defined as the mass in stars that was 
removed by the binary [2]. We find mass deficits of (1.94, 1.56, 1.17) in units of 
the combined black hole mass in the integrations with N = (0.016, 0.065, 0.26) x 
10 6 . These values are of the same order as the mass deficits inferred in giant 
elliptical galaxies [2,21,22]. 

This work was supported by grants AST-0206031, AST-0420920 and AST- 
0437519 from the NSF, grant NNG04GJ48G from NASA, and grant HST-AR- 
09519.01-A from STScI. 
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Figure 3. Initial (dashed line) and final density profiles. 
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